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STATISTICAL ROMBERG EXTRAPOLATION: A NEW VARIANCE 
REDUCTION METHOD AND APPLICATIONS 
TO OPTION PRICING 

By Ahmed Kebaier 

Universite de Marne-La-Vallee 

We study the approximation of E/(Xt) by a Monte Carlo algo- 
rithm, where X is the solution of a stochastic differential equation 
and / is a given function. We introduce a new variance reduction 
method, which can be viewed as a statistical analogue of Romberg 
extrapolation method. Namely, we use two Euler schemes with steps 
5 and 5,0 < (3 < 1. This leads to an algorithm which, for a given 
level of the statistical error, has a complexity significantly lower than 
the complexity of the standard Monte Carlo method. We analyze the 
asymptotic error of this algorithm in the context of general (possibly 
degenerate) diffusions. In order to find the optimal f3 (which turns 
out to be P — 1/2), we establish a central limit type theorem, based 
on a result of Jacod and Protter for the asymptotic distribution of the 
error in the Euler scheme. We test our method on various examples. 
In particular, we adapt it to Asian options. In this setting, we have a 
CLT and, as a by-product, an explicit expansion of the discretization 
error. 

1. Introduction. In many numerical problems — in particular, in math- 
ematical finance — one has to compute, using the Monte Carlo method, 
E,J(Xt), where Xt is a diffusion process. The advantage of using a proba- 
bilistic approach instead of a PDE approach is that one may solve problems 
in high dimension. But on the other hand, the Monte Carlo algorithms are 
much slower and their practical efficiency highly depends on the variance of 
the random variables at hand. This is why variance reduction plays a crucial 
role in the practical implementation. 

There are several classes of methods which are used to reduce variance: 
the control variate approach, the antithetic variate method, moment match- 
ing, importance sampling, conditional Monte Carlo methods. . . . (For more 
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details about variance reduction methods, see [3].) In this paper we intro- 
duce a new variance reduction method which we will call statistical Romberg 
method. This method can be viewed as a control variate method. Roughly 
speaking, the idea is as follows: use many sample paths with a coarse time 
discretization step and few additional sample paths with a fine time dis- 
cretization step. More precisely, in order to construct our control variate, we 
discretize the diffusion (Xt)o<t<T by two Euler schemes with time steps T/n 
and T/m (m <C n). Suppose for a while that M m := E_f(X™) is known. Then 
we construct 

Q = f{X$) - f(Xf) + M m . 

Clearly, we have E/(Xp) = EQ and we wish to simulate Q instead of f{Xj). 
Standard arguments show that for / Lipschitz continuous, we have 

Var(Q) =0(l/m) 

(see Proposition 3.1). So the variance of Q tends to zero as m tends to infinity 
and, consequently, in order to achieve a given accuracy, we need a much 
smaller sample. This significantly reduces the complexity of the algorithm. 
But, on the other hand, we have to compute the quantity M m , and this is 
done again by Monte Carlo sampling. This will increase the complexity of 
the algorithm and we need to find a good balance which guarantees that 
the global complexity decreases. Let N m be the number of Monte Carlo 
simulations used for the evaluation of M m and N n be the number of Monte 
Carlo simulations used for the evaluation of EQ. The question is how to 
choose m, N m and N n . 

In the classical Monte Carlo method, one needs to choose the number n 
of intervals in the discretization and the number N of simulations. The 
parameter n drives the so-called discretization error due to discretization, 
whereas the number N controls the statistical error. For a rational choice 
of N versus n, one may try to minimize the error for a given computational 
time (see [5]) or, equivalently, to minimize the computational time for a 
given total error (see [12]). 

An optimal choice of the number of Monte Carlo samples must be based 
on a precise evaluation of the discretization error. This requires some kind 
of regularity: in [2], the regularizing effect of the diffusion process allows 
to prove that the rate of convergence is 1/n for measurable functions. A 
Hormander type assumption is needed. In the case of general diffusions, the 
same order of convergence can be proved for regular functions only. (See the 
pioneering paper [16] for ^-functions and for ^-functions, see [12]. For 
stochastic differential equations driven by a Levy process, see [8].) 

Our first result shows that, in the context of possibly degenerate dif- 
fusions, the discretization error for ^-functions is at least o(l/^/n) (see 
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Proposition 2.2). For such functions, we give an example for which the dis- 
cretization error is of order l/n a for any a € (1/2,1] (see Proposition 2.3). 

Our second result is a central limit theorem for the statistical Romberg 
method (see Theorem 3.2). This theorem uses the weak convergence of the 
normalized error of the Euler scheme for diffusions proved by Kurtz and 
Protter [11] (and strengthened by Jacod and Protter [9]). 

Based on this central limit theorem, we are able to fix the optimal balance 
between m, N m and N n . It turns out that, for a given error level e = l/n a , 
we obtain m = y/n, N m = n 2a and N n = re 2 " -1 / 2 . With this choice, the com- 
plexity of our algorithm is Csr = C x n 2a+1 / 2 , C > 0, while the complexity 
in the standard Monte Carlo is Cmc = C x n 2a+l . So we have a clear gain. 

Our approach can also be used for the Monte Carlo approximation of 
expectations of functionals of the whole path of a diffusion. In particular, we 
investigate the application to the pricing of Asian options. In this setting, 
the approximation relies on the discretization of the integral of the price 
process over a time interval. It was shown in [17] that the trapezoidal rule 
is one of the most efficient methods for this discretization. We analyze the 
error process, which is of the order 1/ra. We prove a stable functional central 
limit theorem in the spirit of Jacod and Protter [9] . As a consequence of this 
result, we give an expansion of the analytical error, which, in contrast with 
Temam's [17] result, does not use the associated PDE We also use our result 
in order to optimize the choice of the parameters, which are different from 
the ones in the Euler scheme. 

The organization of the paper is the following. We first recall essential 
facts about the Euler scheme, including error evaluations. In Section 3 we 
describe our statistical Romberg method for the Euler scheme. In Section 4 
we apply the idea of statistical Romberg approximation to the discretiza- 
tion of path integrals of a diffusion in the context of Asian option pricing. 
Section 5 is devoted to numerical tests and comparisons. In the last section 
we give conclusions and remarks. 

2. On the discretization error of the Euler scheme. Let (Xt)o<t<T be 
the process with values in M. d , solution to 

(1) dX t = b(X t )dt + a(X t )dW t , X = xeR d , 

where W = (W 1 , . . . , W q ) is a q-dimensional Brownian motion on some 
given filtered probability space B = (O, (Ft)t>o, P)- (-^i)t>o is the stan- 
dard Brownian filtration. The functions b : M rf — ► R rf and a : M. d — ► M. dxq are 
continuously differentiable and satisfy 3 Ct > 0; Vx, y G M. d , we have 

\b(x) - b(y)\ + \a{x) - a(y)\ < C T \y - x\. 

We consider the Euler continuous approximation X n with step 5 = T/n 
given by 

dX? = b(X Vn(t) )dt + a(X nn{t) )dW t , Vn (t) = [t/6}6. 
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It is well known that the Euler scheme satisfies the following properties 
(see, e.g., [6]): 



(2) Vp>l 

(3) Vp>l 



Esnp \X t -X?f 
te[o,T] nPl 1 



K p (T) > 0. 



E sup \X t \ p + E sup \X?\ P < K'JT), K'JT) > 0. 
te[o,T] te[o,T] 



2.1. Stable convergence and the Euler scheme error. We first recall basic 
facts about stable convergence. In the following we adopt the notation of 
Jacod and Protter [9] . Let X n be a sequence of random variables with values 
in a Polish space E, all defined on the same probability space (fi, J-, P). Let 
(p,,J-,P) be an extension of (Cl,J-,P), and let X be an E- valued random 
variable on the extension. We say that (X n ) converges in law to X stably 

and write X n s ^$ r x, if 

E (Uh(X n )) — ► E (Uh(X)), 

for all /i : E — > R bounded continuous and all bounded random variable £/ 
on (fi, J 7 ). This convergence, introduced by Renyi [15] and studied by Aldous 
and Eaglson [1], is obviously stronger than convergence in law. The following 
lemma will be crucial: 

If V is another variable with values in another Polish space F, we have 
the following result: 

Lemma 2.1. Let V n and V be defined on with values in another 

metric space E' . 



IfV n ^V,X n s ^X, 



then (V n ,X n ) s ^ y (V,X). 



This result remains valid when V„ = V . 



For a proof of this lemma, see [9]. 

Note that all this applies when X n , X are R d - valued cadlag processes, 
with E = B([0,T],R d ). Now assume that 



<p(Xt 



fbx{X t ) <Tu(X t ) 
b 2 {X t ) a 21 {X t ) 



<T 2q {X t ) 

o- dq (X t )J 



\b d (X t ) a dl (X t ) . 
then the SDE (1) becomes 
(2) dX t = (p(X t )dY t . 



I dt \ 
dW} 



and dYt 



\dW t q J 
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The Euler continuous approximation X n with step 5 = T/n is given by 

dX? = <p(X Vn{t) )dY t , rj n (t) = [t/8]5. 

The following result proven by Jacod and Protter [9] is an improvement on 
the result given by Kurtz and Protter [11]. 



Theorem 2.1. With the above notation we have 
JnU n =: ^{X n - X) S ^ U, 
with U a d- dimensional process satisfying 



q+l d 

(3) dUt = J2T,Vk(Xt 
j=ik=i 



9+1 



U t k dY{-J2<P kl (Xt) dN l t 



E/g = 



(<f^ 3 is the partial derivative of ip % i with respect to the kth coordinate), and 
N is given by 



N u = 0, 
N jl = 0, 



l < « < g + 1, 
i<j<g + i, 

N ij = ^, 2<i,j<q + l, 

where (-B y )i<jj<g is a standard (q) 2 -dimensional Brownian motion defined 

on an extension $= (^l,^,(T t )t>o,P) of the space (0, T, (^Pt)t>o, P)> which 
is independent ofW. 



We will need the following property of the process U . 

Proposition 2.1. Under the assumptions of the above theorem, we have 

M(U T /F T ) = 0. 



Proof. Consider the unique solution of the d-dimensional linear equa- 
tion 



(4) 



9+1 fT 



where (p' 1 is a d x d matrix with ((f'^ik = ip'u ' ■ From Theorem 56, page 271, 
in [14], it follows that 



(5) 



<?+l .t 

Ut = _J2S t £rWx t )^X t )dNl 

.7=1 J ° 
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If Z is a bounded J^-measurable r.v., we have 

E (U T Z) = " E E [ Z£ t J o £ t V {X t )ip{X t )dNij = 0, 
as can be seen by representing Z£t as a stochastic integral w.r.t. W. □ 

2.2. TTie discretization error. In the following we focus on the discretiza- 
tion error given by the bias 

(6) e n :=Ef(X2)-Ef(X T ), 

where / is a given function. Talay and Tubaro [16] prove that if / is suffi- 
ciently smooth, then e n ~ c/n with c a given constant. A similar result was 
proven by Kurtz and Protter [12] for a function / £ c to^. The same result 
was extended in [2] for a measurable function /, but with a nondegeneracy 
condition of Hormander type on the diffusion. In the context of possibly de- 
generate diffusions, the discretization error for functions which are not ^ 3 
is not yet completely understood. For a Lipschitz-continuous function /, the 
estimate |e n | < ^= follows easily from (2). The following proposition and 

the example below focus on the rate of convergence of e n for ff 1 functions. 

Proposition 2.2. Let f be an R d -valued function satisfying 

(7) \f(x)-f{y)\ <C{l + \x\ v + \y\ p )\x-y\ for some C,p > 0. 

Assume that F(X T $ V f ) = 0, where V f := {x £ R d \f is differ entiable at x}, 
then 

lim \fne n = 0. 

n— >oo 

Proof. We have, with probability 1, 

vW(*?) - f{X T )) = v^V/(X T ) • U$ + R n , 

with 

R n = ^n~\U2\e(X T ,U£) and e(X T ,U$)$0. 

p 

It follows that R n — ► 0, since (v^l^rl) ^ s tight- Consequently, we deduce, 
using (13), (2), (3) and Theorem 2.1, that 

^e n ^M(Vf(X T )-U T ), 
and using Proposition 2.1, it follows that 

i(v/(x T ) • u T ) = o, 

which completes the proof. □ 
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The following example proves that, for a G (1/2,1], there exists a < ?f 1 
function with bounded derivatives and a diffusion X such that 

(7) n a e n ^C f (T,a), 

where Cf(T,a) is positive. In other words, the rate of convergence can be 
l/n a for all values of a G (1/2, 1]. 

Example 2.1. Consider the bi-dimensional diffusion Z= (X, Y) satis- 
fying the following SDE: 

dX t = -X t /2dt-Y t dW t , 

(8) 

dY t = -Y t /2dt + X t dW t 

and the map f a :z = (x,y) \-> \\z\ 2 — l\ 2a . The solution of (8), subject to 
Zq = (cos 9, sin#), is given by Z t = (cos(0 + Wt),sm(9 + Wt)). We assume 
that 9 G [0,2-71"], so the diffusion Z lives on the unit circle. 

Proposition 2.3. Lei Z n be the Euler scheme associated with Z . For 
a G [1/2, 1], we have 

(9) n Q E(/ a (Zf) - / a (^)) - (2t) a E|G| 2a , i > 0, 

where G is a standard normal r.v. 

Proof. We have, since f a vanishes on the unit circle, 
n a E(f a (Z?) - f a {Z t )) = n a E\\Z?\ 2 - l\ 2a 

= n a E\[Z t + Z?]-[Z?-Z t ]\ 2a . 

Using Theorem 2.1, 

>t\ , 

where U = (£/ , U 2 ) is given by 



(10) n a E(f a (Z?) - f a (Z t )) - 2 2a E|Z t • C/ 4 | 2 ° 



dE/* = -\fj\dt- U? dW t + -j^X t dB t 



(11) 

dC/ 2 = - if/ 2 dt + E# dW t + ^Y t dB t 

and is a standard Brownian motion independent of W. The solution of (11) 
is given by 

(12) 



which completes the proof. □ 
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3. The statistical Romberg method. Before introducing our algorithm, 
we recall some essential facts about the Monte Carlo method. In many ap- 
plications (in particular, for the pricing of financial securities), the effective 
computation of Kf(Xx) is crucial (see, e.g. [13]). The Monte Carlo method 
consists of the following steps: 

• Approximate the process (X t )o<t<T by the Euler scheme (X t n )o<t<r 5 with 
step T/n, which can be simulated. 

• Evaluate the expectation on the approximating process f(Xj) by the 
Monte Carlo method. 

In order to evaluate E/(X^) by the Monte Carlo method, N independent 
copies /(Xp j)i<i<jv of f{Xj,) are sampled and the expectation is approxi- 
mated by the following quantity: 

1 N 

4 £/(*?,*)• 

i=X 

The approximation is affected by two types of errors. The discretization error 
e n , studied in the above section, and the statistical error f n,N — E/(Xy), 
controlled by the central limit theorem and which is of order 1/y/N. An 
interesting problem (studied in [5] and [12]) is to find N as a function of n 
so that both errors are of the same order. 

The following result highlights the behavior of the global error in the 
classical Monte Carlo method. It can be proved in the same way as the limit 
theorem given in [5]. 

Theorem 3.1. Let f be an M. d -valued function satisfying 

(13) \f{x) -f(y)\ <C(1 + \x\ p + \y\ p )\x-y\ for some C,p > 0. 

Assume that ¥(Xx ^ T>A = 0, where T> y := {x G is differentiable at x}, 
and that for some a G [1/2, 1], we have 

(14) lim n a e n = C f (T,a). 

n— »oo 

Then 

n ° E f( X h) - E K x t)) ^"JG + C f (T, a), 
with a 2 =Var(/(Xr)) and G a standard normal. 

A functional version of this theorem, with a = 1, was proven by Kurtz 
and Protter [12] for a function / of class ^ . We can interpret the theorem 
as follows. For a total error of order l/n a , the minimal computation effort 
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necessary to run the Monte Carlo algorithm is obtained for N = n . This 
leads to an optimal time complexity of the algorithm given by 



Recall that the time complexity of an algorithm A is proportional to the 
maximum number of basic computations performed by A. 

3.1. The Euler scheme and the statistical Romberg method. It is well 
known that the rate of convergence in the Monte Carlo method depends on 
the variance of f{XJf), where Xj, is the Euler scheme of step T/n. This is a 
crucial point in the practical implementation. A large number of reduction 
of variance methods are used in practice. Our algorithm proposes a control 
variate reduction of variance. Its specificity is that the control variate is 
constructed itself using the Monte Carlo method, applied to the same dis- 
cretization scheme, but with a step m which is specifically lower than the 
approximation step n (using two discretization steps is an idea which al- 
ready appears in Romberg's method). Let us be more precise. We fkm«n 
and we denote 



where M m = E/(Xy) and we suppose for a while that M m is known. Note 



so that f(X™) — M m appears as a control variate. 

Consider a function / : R d — > M. d which is Lipschitz continuous of constant 



(13) C M c = Cx (nN) = Cxn 



2a+l 



with C some positive constant. 



Q = f(X$) - f{XJp) + M, 



that 



E(Q)=E/(X£) 



[/] lip , that is, [/] Up = sup u _^ IJ ^_^[ 



Proposition 3.1. Under the above assumptions, we have 



(14) 



0", 



.2 

Q 



Var(Q) = 0(l/m). 



Proof. We have 



a Q = \\Q-EQ\\ 2 

<\\f(x®-f{xF)h 



<[/]ii P sup \\X t -X?\\ 2 + sup \\X t -X™\\ 2 . 



Ue[o,T] te[o,T] 
Using (2), we deduce that 3K' > such that 
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which completes the proof. □ 

Inequality (14) shows that the variance of Q is significantly smaller than 
the variance of f(X^), so that Q appears as a good candidate for the re- 
duction of variance method. However, computing EQ supposes to compute 
M m = E/(Xj?) in the first place, and this is also done by the Monte Carlo 
method. So there is a certain extra quantity of computation to be done. In 
practice, the sample paths used for the computation of M m will be indepen- 
dent of those used for the computation of Q. 

In the following we make a complexity analysis which permits to choose m 
as a function of n in order to minimize the complexity of the algorithm which 
leads to EQ. We will prove that, with such a choice of m, the complexity 
of the algorithm for EQ is significantly smaller than the complexity of the 
standard Monte Carlo method for ~Kf{Xj). 

Let us present the algorithm for EQ. 

3.2. Central limit theorem. In the following we assume that the param- 
eters of the statistical Romberg method depend only on n. That is, 

m = nP, /3e(0,l), N m = n<\ 71 > 1, N n = n^, 72 > 1. 

We can now state the analogue of Theorem 3.1 in our setting. The statistical 
Romberg method approximates E/(Xy) by 

i=l i=l 

where Xjf is a second Euler scheme with step T/nP and such that the Brow- 
nian paths used for X^, and Xjf have to be independent of the Brownian 
paths used in order to simulate Xlf . Here the quantity p-ES/^T^) 
must be viewed as an approximation for M m . 

Theorem 3.2. Let f be an W 1 -valued function satisfying 

(15) \f(x)-f(y)\ <C(l + \x\ p + \y\ p )\x-y\ for some C,p > 0. 

Assume that F(X T V f ) = 0, where V f := {x £ R d ; f is differentiable at x}, 
and that for some a € [1/2, 1], we have 

(16) lim n a e n = C f (T,a). 

Then, for 71 = 2a and 72 = 2a — (3, we have 

n a (V n -Ef(X T ))^a 2 G + C f (T,a), 
with a\ = Var(/(-Xy)) + V "ar(V f '(Xt)Ut) and G a standard normal. 
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Lemma 3.1. Under the assumptions of Theorem 3.2, for all 7 > 0, 

( 15 ) - f{x T ,)-nf{xjf) - f(x T )) =^(0,(7?), 

(16) a? = Var(V/(X T ) • t/ T ) 
and U the process on B given by (3). 

Proof. If we set 



Zf = f(xf) - f{X T ) - E(f(xf) - f(X T )), 



then we have 



E 



1 + — ( — n% |Z^| 2 + nmC n {Lu) 
n 7 V 2 



where 



Property (2) ensures the existence of a constant K3 > such that 

, \ 1 K 3 u 3 

1^)1 <g^. 

We have, with probability 1, 

n^ 2 (f(xf) - /(X r )) = n^ 2 V/(X T ) • < + i^, 

with 

R n = n^ 2 \uf\e(X T ,uf) and e(X T ,utf) $0. 

From the tightness of (n^^ 2 \U^ \) n , it follows that R n 0, then, according 
to Lemma 2.1 and to Theorem 2.1, 

(17) n^ 2 (f(xf) - f(X T )) => V/(X T ) • U T . 

Using (13), it follows from property (2) that 

Ve > mpE\n p/2 (f(X^) - f{X T ))\ 2+£ < 00. 

n 

Since F(X T £ V f ) = 0, we deduce, using (17), that 
E(n^ 2 (f(xf) - f(X T ))) k ^E(Vf(X T ) ■ U T ) k < 00 with k G {1,2}. 
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Consequently, 

(18) n^E\zf\ 2 — >Var(V f(X T ) • U T ) <oo. 

Since 7 > 0, we see that n 7 EC„(w) — ► and we conclude that 



E 



cxp 



n (( 7 -/3)/2) 



cxp 



-It 



Var(V/(X T ) • 17 T ) 



which completes the proof. □ 



Lemma 3.2. Under the assumptions of Theorem 3.2, for all 7 > 0, 
1 ™ 7 

£ /(X£, fc ) - f(X$ k ) - E(/(X£) - /(X/)) =^#(0,^), 



((7-/3)/2) 



7? 



fe=l 



where a\ = Var(V/(X T ) • U T ). 

PROOF. We have 

1 nl 

n (( 7 -fl/2) E /(^t,.) " " E(/(JC?) - /(X^)) 

k=l 

1 _ 1 



n (( 7 -/3)/2) "'J'.* 
fc=l 



E ^?,A 



n ((7-/3)/2) 



k=l 



By (18), it follows that 



(19) 



E 



1 



n ((7-/9)/2) 



fc=l 



n% [Z£] 2 



0. 



The announced result follows from the above lemma. □ 

Proof of Theorem 3.2. For 71 = 2a,72 = 2a — (3, we have 
n a (V n -E/(X T )) = + V* + 

where 

1 0, 



(20) # = — -E/(*£ „ 

n i=i 

(21) V^—Y^ f(X^) - f(X^) - E(/(X?) - f(xf)), 



(22) y n 3 = n Q (E/(^)-E/(X T )). 



i=i 
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Properties (2) and (3) guarantee that the Lindeberg-Feller theorem applies 
here (same argument as in [5]). That is, 

On account of Lemma 3.2, it is obvious that 

V 2 =► AT(0, Var(V/(A T ) • U T )). 
Finally, by using the assumption (14), we complete the proof. □ 

3.3. Complexity analysis. As in the Monte Carlo case, we can interpret 
Theorem 3.2 as follows. For a total error of order l/n a , the minimal com- 
putational effort necessary to run the statistical Romberg algorithm applied 
to the Euler scheme (with step numbers n and m = n^) is obtained for 

(23) N m = n 2 a and N n = n 2a ~ p . 

Since the only constraint on (3 is that j3 S (0, 1), we will choose the optimal (3* 
minimizing the complexity of the statistical Romberg algorithm. The time 
complexity in the statistical Romberg method is given by 

C SR = C x (mN m + (n + m)N n ) with C > 

= C x ( n P +2a + {n + nP)n 2a -P). 

Simple calculations show that (3* = 1/2 is the optimal choice which mini- 
mizes the time complexity. 

So the optimal parameters in this case are the following: 

m = n 1/2 , N m =n 2a and N n = n 2a ~ 1/2 , 

and the optimal complexity of the statistical Romberg method is given by 

C SR ^C7xn 2Q+1 / 2 . 

However, for the same error of order l/n a , we have shown that the optimal 
complexity of a Monte Carlo method was given by 

C MC = Cxn 2a+1 , 

which is clearly larger than Csr. So we deduce that the statistical Romberg 
method is more efficient. 



4. Statistical Romberg method and Asian options. The payoff of an 
Asian option is related to the integral of the asset price process. Comput- 
ing the price of an Asian option requires the discretization of the integral. 
The purpose of this section is to apply statistical Romberg extrapolation 
to the approximation of the integral and to carry on a complexity analysis 
in this context. This will lead us to prove a central limit theorem for the 
discretization error, which can be viewed as the analogue of Theorem 2.1 
(see Theorem 4.1). 
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4.1. Trapezoidal scheme. Let S be the process on the stochastic basis 
B = (fi,.F, (Ft)t>o,P) satisfying 

7 C 

(24) — - = rdt + CTdW t with t€ [0, T], T > 0, 

St 

where tr and r are real constants, with a > and (W^g^T] is a standard 
Brownian motion on B. The solution of the last equation is given by 

(25) S t = S exJ(r-^\t + aW t 



1 f T 
It = ^ I S u du. 
o 



We set 

Trr = 

T 

Let / be a given real valued function. Our aim will be to evaluate 

U(S,T) = e- rT Ef(S T ,I T )- 

In a financial setting, if /(x, y) = {y — K) + , H(S, T) is the price of an Asian 
call option with fixed strike K. In this case there is no explicit formula 
that gives the real price. So, the computation of this price, by a probabilis- 
tic method, requires a discretization of the integral iy. There are several 
approximation schemes used in practice. One of the most efficient is the 
trapezoidal scheme defined by 

( 26 ) I T = j,2^S tk _ 1 [l + — + a j, 

where 5 = ^ and tk = ™ = 5k. We call it trapezoidal because it is closely 
related to the trapezoidal approximation of the integral 

Note that St k has an explicit expression so we can simulate it without dis- 
cretizing the SDE. The following result is proved in [17]. 

Proposition 4.1. With the above notation, there exists a nondecreasing 
map K(T) such that, V p > 0, 

V(2p) K(T) 
~ n 



E( sup |/ t n -/t| 2p 
te[o,T] 
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4.2. Stable convergence of the trapezoidal scheme error. In the following 
we prove a functional CLT analogous to Jacod and Protter's theorem (see 
Theorem 2.1 above). 



Theorem 4.1. Let J t = ^ Jq S u du,t £ [0, T], and J n be the trapezoidal 
discretization associated with J: 



[t/5] 



S ^ „ A r5 W tk - Wt h _, 



T 

k=i 



We have 



(27) n(J - J n ) s ^l v X , 

where x is the process defined by 

where B' is a standard Brownian motion on an extension B of B, which is 
independent of W. 

We have the following elementary lemma. 

Lemma 4.1. If H is deterministic satisfying Jq ds < oo, then we have 



f lM Hsds _^^ r Hsds 

Jo o 2 Jo 



and 

with t s (s) = t A ([s/5}5 + 5) -s. 

Proof. We sketch the proof for completeness. By a density argument, 
we may assume that H is piecewise constant: H s = Ci for Tj_i < s < Ti, 
where = To < • • • < = T and (q) are constants for i = 0, . . . , k. It follows 
that 

f T -f H "» = £jC' r. W * - £ f Pi*. " V) = [ f *, 
since it is easy to check that 

l f v { s , y-x 

- / T s (s)as — > — - — as n — > oo. 
o ' 



.r 
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The second assertion is obtained in the same way, but by using that 

1 [V ( tAs) e \ 2 , y-x 

2 — o as — > as n — > oo, 



which completes the proof. □ 

Proof of Theorem 4.1. We have 

n ? (5 m 5 2 r [t/5] 

n(Jt ~ Jt) = y J o S u du-n\^-Y, St fc _! + yfYl S ^-i 

5a [t/S] \ 
(28) +^ y £S th _ 1 (Wt k -Wt k _ 1 )j 

+ y s ws\s (i + \(t - + a -{w t - w ms )j . 



It follows that 



(29) n(J t - J?) =Al- T -f S [u/S]5 du-^J* S [u/S]s dW u , 



with 



1 rt 



At = 6 Jo ( Su ~ S [ u / S ] 5 "> du - 



Note that, by using (24), we obtain 



with 



and 



Si r rt ru 
A t ' = - I I Ssds du 



S Jo J[u/S]S 

A f 2 = T f f U S s dW s du, 
Jo J\u/8]8 



and we have 



s 1 r 1^-5 /"** f u r /"* f u 

A t ' = -y I I S s dsdu + - / S s dsdu, 

° t^i Jt k-i Jtk-i 6 J[t/S\8 J[t/S]8 

- / (tA([s/5]5 + 5)-s)S s ds, 



o 



/ T s (s)S s ds. 
Jo 
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In the same manner we can see that 

I*/ 5 ] „t k fU g- ft fU 

4' 2 = tE/ / S s dW s du + - / S s dW s du. 

j£ A fe -i /t fc _i « ^[t/<5]5 J[t/S]S 

The integration by parts formula gives 

ft* /•« /-tfc /-tfc 



/ / S s dW s du=(t k -t h _ 1 ) S s dW s - (s-t k )S s dW s 

f tk ([s/S]S + S-s)S s dW s 
Jtk-i 



and 



S s dW s du = (t- [t/S] 5) / 5 S dW s 

[t/S]S J[t/8}5 J[t/5}5 

( S -[t/5]S)S s dW s 

WW 



[ (t - s)S s dW s 
J\tlS\S 



'WW 
We deduce that 

CT [t/5] /*k (J ft 

A t' 2 = lE / ([ s /<5]5 + 6 - s)S s dW s + - (t- s)S s dW s , 

o /t t _i o J\tlS\5 



k=l 



G - [\tA([s/S]S + S)-s)S s dW s 
o Jo 



^ f T s {s)S s dW s . 
o Jo 



It follows that 



n(J t - J") = - f r s (s)S s ds + -[ l Tg (s)S s dW s 
o Jo o Jo 



a 



S\s/S]S ds - — / Su/s]s dW s 
i Jo £ Jo 

We deduce that 

n(J t -jn=Bf + X 5 t +C s t , 

with 



>5 r 



B °t = 2 I- (2^-1)^^, 
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- S[ s/S ]s)dW s . 



According to the above lemma, we obtain that sup^rQ ^i Bf — > a.s. It is ob- 
vious that sup tg [ 0j r] W - ► 0. The only point remaining concerns the behavior 
of x S - 

In virtue of Theorem 2.1 of [7], if we prove that, for all t S [0, T], we have 

2 t 

12 Jo 

then the process x & will converge stably in law to the process x °f (27). 
According to the above lemma, we have 

(x 5 ,X S )t = £ f o (2 T -M - l) 2 Si dW s — g jf Si ds a.s. 

and 

(x\ W) t = \ J* (2^1 -iy s ds^0 a.s., 

which completes the proof. □ 

4.3. Statistical Romberg method and CLT. In order to evaluate e~ rT Kf(ST, 
It), we use the idea of statistical Romberg approximation: 

• compute an approximation E\ of e~ rT Kf(ST,lT ) by a Monte Carlo 
method 

e -rT " 71 

E n = Yl f( S T,h Ir,i), 

i=l 

• compute 

i=l 

Recall that the samples used, in order to construct 

and ((^i,!^),^,!^)), 

are independent. The question now is how to choose /3,7i and 72. Admit- 
tedly, we can choose the optimal parameters given in the Euler scheme case, 
but the following result proves that, in the specific case of the trapezoidal 
scheme, the optimal parameters are different. 



(30) 
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Theorem 4.2. Let f be an R 2 -valued function satisfying 

\f(x, yi ) - f(x,y 2 )\ < C(l + \x\p + \ Vi \p + \y 2 \P)\yi - y 2 \ 

for some C,p>0. 
Assume that P((S T ,I T ) i %) = 0, with 

@ f ■= {(x,y) e R + x R + | d 2 f(x,y) exists}, 

where d 2 f denotes the partial derivative of f w.r.t. y. Then for all (3 G (0, 1), 
if 7i = 2 and 72 = 2 — 2/3, w;e aaue 

n(£7„ - E/(5 r , It)) =>> <r 2 G + E(5 2 /(5 T; I T )xr), 

where o\ = Var(/(Sr, iy)) + V&r(d 2 f(ST, It)Xt)> X ^ the limit process on 
B given in Theorem 4.1, and G a standard normal. 

Remark 4.1. The assumptions on / in the above theorem are satisfied 
in the case of typical Asian options: 

f(x,y) = (y- K) + , f(x,y) = (K -y)+, f(x,y) = (y-x)+. 
Lemma 4.2. Under the assumptions of Theorem 4.2, for all 7 > 0, 
1 nl n 

X] f(ST,kjT,k) - f(ST,k,lT,k) 



n ((7/2)-/3) 



fc=l 



-E(f(S T ,I T ) - f(S T jf)) ^A/(0,<r 2 ), 
where a 2 = Var(d 2 f(S T ,I T )x T )- 

PROOF. If we set 

Hf := f(S T ,I T ) - f(S T ,lf ) - E(f(S T ,I T ) - f(S T ,lf)), 



then we have 



E 



vu 



cxp 



„(( 7 /2)-/3) 



(30) 



k=l 
-U 2 



1 + -7 hr VarK(/(5 T , It) - /(5 T ,i? ))] 



+ n 7 EC^(a;) 



with 
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Proposition 4.1 ensures the existence of a constant K(T) > such that 

We have, with probability 1, 

n^f(S T ,lf) - f(S T ,I T )) = n^d 2 f(S T ,I T )(lf - I T ) + R n , 

with 

R n = n p \lf -ItHSt, It, if) and £(S T ,I T ,lf)^0- 

From the tightness of {nP\lf — Ir|)n> it follows that R n ^> 0. Consequently, 
according to Lemma 2.1 and to Theorem 4.1, we obtain that 

(31) n P(f(ST,lf) - f(S T ,I T ))^d 2 f(S T ,I T )x T . 

With our assumption (30) on /, it follows from Proposition 4.1 that 



supE|n /3 (/(5 T ,/T) - /(StJt))^ < °° with e > 0, 



rA\|2+e 



so we obtain 



(32) 



E(nP(f(S T ,I T )-f(S T ,lf))) k 



E(d 2 f(S T ,lT)x T ) <°o V0<A;<2. 



Hence, we deduce that 



E 



exp 



(- 



m 



I n (( 7 /2)-/3) " 
\ fe=l 

which completes the proof. □ 



exp 



-u 



Y & i(d 2 f(S T ,I T )x, 



Remark 4.2. It follows from the proof of the above lemma that 
(33) lim nM(f(S T ,m - f(S T ,I T )) ^E(d 2 f(S T , I t ) X t)- 

n—>oo 

This gives us an expansion of the discretization error in our setting. Note 
that similar expansions are given in [17] for less regular functions. The 
advantage of our approach is that we do not need the associate PDE, so 
that our expansion is more explicit. 



The proof of the following result is a consequence of the above lemma. 
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Lemma 4.3. Under the assumptions of Theorem 4.2 and for all 7 > 0, 
we have 

1 nl 

„(( 7 /2)-/3) 5^ /(^r.fc, Er,k) - f( S T,k, Ir,k) 
fc=i 

- E(/(5 T , i?) - /(5r, ^)) AT(0, 
wi/i of = Yar(d2 /(St, It)x t )- 

Using Lemma 4.3, Theorem 4.2 can be proved in much the same way as 
Theorem 3.2. Equality (33) will be used instead of the assumption (14) on 
\\m n ^ 00 n a e n in Theorem 3.2. 

4.4. Complexity analysis. As in the Euler scheme case, one can interpret 
the above theorem in the following way. For a total error of order 1/n, 
the minimal computational effort necessary to run the statistical Romberg 
method applied to the trapezoidal scheme (with step numbers n and m = n^) 
is obtained for 

(34) N m = n 2 and N n = n 2 ~ 2fi . 

Since the only restriction on j3 is that (3 E (0, 1), we will choose the optimal 
/3* minimizing the complexity of the statistical Romberg algorithm. In this 
case the time complexity in the statistical Romberg method is given by 

C S r = C x (mN m + (n + m)N n ) with C > 0, 

= Cx(n' 3+2 + (n + n' 3 )n 2 - 213 ). 

Simple calculations show that (3* ~ 1/3 is the optimal choice which mini- 
mizes the time complexity. 

So the optimal parameters in this case are 

m = n l/3 ? N m = 2 and N n = n 4/3 , 

and the optimal complexity of the statistical Romberg method is given by 

CsR~Cxn 7 / 3 . 

But according to Proposition 4.1, Theorem 3.1 remains valid if we change the 
Euler scheme by the trapezoidal one. Hence, for the same error of order 1/n, 
the optimal complexity of a Monte Carlo method applied to the trapezoidal 
scheme with step number n is given by 

Cmc = C x n 3 , 

which is clearly larger than Csr. So we deduce that the statistical Romberg 
method is more efficient. 
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5. Numerical tests and results. We test the efficiency of the statistical 
Romberg method to reduce the time complexity for the degenerate two- 
dimensional diffusion given in the example of Section 2 (tests concerning 
Asian options are given in [10]). 

Consider the bi-dimensional diffusion 

Z t = (cQB(0 + Wt),sm(P + Wt)), 0€[O,27r], 

and the map 

f a :z = (x,y)^\\z\ 2 -l\ 2a , ae [1/2,1], 
given in the example of Section 2. Note that in this case 

Var(/ a (2r)) = 0, 

and this is because the bi-dimensional diffusion Z = (X,Y) lives on the 
unit circle. Consequently, in order to obtain the optimal parameters given 
by Theorems 3.1 and 3.2, we consider g a (x,y) = f a (x,y) + x, instead of 
f a (x,y). This choice leads us to a nonvanishing variance of g a (Zx) and to a 
discretization error which is of order l/n a , a E [1/2,1]. In the following we 
set 

• MC method: the algorithm using a Monte Carlo method to approximate 
E(g a (Z T )) by 

1 N 

i=l 

• SR method: the algorithm using a statistical Romberg method to approx- 
imate E(g a (Z T )) by 

i N n 1 N m 

( 35 ) *r£M3?,i) -9a(z^)} + — Y,g a (z^). 

To compare both methods, we use the methodology proposed by Broadie and 
Detemple [4]. Their idea is that, for a given set of parameters of the con- 
cerned diffusion, one of both algorithms will give better results. So they 
propose to test the algorithm on a large set of parameters chosen ran- 
domly. Proceeding along this line, we produce randomly M = 200 values 
for Zq = (Xq,Yq). Then, for each method, we compute the speed and an 
error measure. Speed is measured by the number of simulated values com- 
puted per second (the computations were done on a PC with a 2.00 GHz 
Pentium 4 processor). The error measure is given by the root-mean-squared 
error, which is defined by 



(36) RMS 



1 M 

— ^(Real value — Simulated value) 2 , 

i=l 
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and the real value is given by the formula K(g a (ZT)) = cos(# — T/2). 

Our algorithm proceeds as follows. We fix the number of steps, say, n = 80, 
in the Euler scheme. We compute by the Monte Carlo method (resp., the 
statistical Romberg method) the 200 simulated values. Then, we produce 

according to (36) RMS (MC ' n=200) and we compute Speed (MC ' n=200) , so we 
have a couple of points 

(RMS (MC,n=200) ) Speed (MC,n=200) ) 

This point is plotted on Figure 1. So, each given n produces a couple of 
points: 

(MC (MC,n) )Speed (MC,n) )/rhis 

gives the continuous curve in Figure 1 . 
The line marked by squares is produced in the same way, using the statistical 
Romberg method. 

Tests are done for a = 1/2. Note that, although g a is not c to l in that case, 
Theorem 3.2 can be extended to this specific example, as we can handle the 

difference ffi/ 2 (^T ) ~9i/2{ z t)- 

Let us now interpret the curves. The fact that the SR curve is higher 
than the MC one means that, given an error e, the number of values com- 
puted in one second (with this error) by the statistical Romberg algorithm 
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Table 1 




Time 


complexity reduction for 


a = 1/2 


RMS relative error MC method speed 


SR method speed 


HT 1 


45.405 


102.718 


9-l(T 2 


23.599 


85.876 


8- 10~ 2 


18.647 


70.712 


7- 1CT 2 


9.219 


49.179 


6- 1CT 2 


5.883 


29.234 



is larger than the number of values computed by the Monte Carlo method. 
Note anyway that, for a large s (which corresponds to a small number of 
steps n), the differences between the two methods is less important. But as e 
becomes small (n becomes large), the difference becomes more significant. 
In Table 1 we compare the speed of the Monte Carlo method and the speed 
of the statistical Romberg one, for a fixed RMS-error. We note that, by 
using the statistical Romberg method and for an RMS-error fixed at 10 , 
one increases the speed by a factor of 2.26. For a small RMS-error fixed 
at 6 • 10 -2 , the speed gain reaches a factor of 4.96. 

6. Conclusion. The statistical Romberg algorithm is a method that can 
be used in a general framework: as soon as we use a discretization scheme 
for the diffusion (X t )o<t<T in order to compute quantities such as Kf(Xx), 
we can implement the statistical Romberg algorithm. And this is worth it 
because it is more efficient than a classic Monte Carlo method. 

In financial applications, it is sometimes essential to be able to price a 
given product on the market as soon as possible and this by setting a margin 
of error that one can tolerate. In this case the statistical Romberg method 
equipped with its parameters allowing complexity reduction is faster than 
the standard Monte Carlo one. 
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